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1.0 INTRODUCTION 

Four programs have been developed to perform the calcula- 
tions of flow field in the nose region of a blunt body at 
arbitrary speed and altitude. These programs differ from 
each other in their ability to consider either thin shock 
or thick shock conditions and in the use of either ideal, 
equilibrium air or nonequilibrium air chemistry. No attempt 
has been made to incorporate all four programs into a single 
one, and to use option cards to choose the particular program 
for use. This effort would be quite time-consuming because 
in each of the four programs there exists a number of options, 
viz., the choice of three orthogonal coordinate systems, four 
different body surface conditions in temperature and cata- 
lycity, and simple or multicomponent transport properties. 

It is felt that in order to achieve better computational 
efficiency and lower computer storage requirement, different 
programs should be used for different purposes. For clarity 
of further discussion, these four programs are termed as 

I. Non-reacting thin-shock code 

II. Non-reacting thick-shock code 

III . Reacting thin-shock code 

IV. Reacting thick-shock code 

These programs have been completely checked out and applied 
successfully to the flow field computations described in Part 
1 of this report. This part of the report will concentrate 
on the use and restrictions of these programs. In general, 
accurate solutions can be obtained for a blunt body flow 
with or without finite-rate chemical reactions for an alti- 
tude less than 400 Kft and a speed less than 26 Kft/sec, 
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or for a free stream Reynolds number, Re^/Rj^ , larger than 

one. However, it should be mentioned that since the computer 
execution time for Programs III and IV is very high, in 
the range of 3 - 10 hours on a UNIVAC 1108, present experience 
with these programs is gained from several selective cases. 
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2.0 PROGRAM DESCRIPTION 

These programs are written primarily for solving Eqn (10) 
of Part 1 by means of a finite-difference method. The 
method of solution is described in Section 4, of Part 1. 

A description of the non-reacting thin and thick-shock 
codes is first given which is followed by the description 
of the reacting thin and thick-shock codes. The difference 
between the programs of thin and thick codes is indicated 
when it occurs. 

MAIN subroutine is used to read in the parameters of the 
problem, i.e., the number of mesh points and time steps, the 
free stream conditions , etc. The subroutine BODY defines 
the body geometry which is presently limited to analytic 
equations. The initialization of the flow field computation 
is made in subroutine INIT. INIT is not called in the thick- 
shock code for which the necessary initial flow field is 
obtained in subroutine CONVR1, along with the conversion 
of physical flow properties into the vector U of Eqn (14) 
and (15) . The metric coefficients of the coordinates used 
are determined in subroutine REGION. The time-marching 
forward procedure of Eqn (23) is programmed in subroutine DIFF. 
The boundary conditions as specified in equations in Section 
2.4 are imposed to the computed results by subroutines 
SHOCK and WALL. The subroutine SHOCK is not used in the 
thick-shock code since the outer boundary has the free 
stream conditions. The constitutive equations defining 
the stress tensor and heat flux vector are coded in sub- 
routines STRESS and HEAT, respectively. Thermodynamic 
properties of the flow are obtained by calling the subroutine 
THERM. The input to THERM is obtained from subroutine 
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C0NVR2 , which converts the flow vector U to the physical 
properties. Finally, the subroutine OUTPUT prints out 
the results. 

The program structure of the reacting thin and thick-shock 
codes is similar to that of non-reacting codes. Therefore, 
only the additional subroutines are given. These subroutines 
are needed to read information about the thermodynamic co- 
efficients and transport properties of the mixture, and the 
chemical reactions. Subroutines INPUTT and INPUTR are 
responsible for this purpose. The initial flow field can be 
defined with a better accuracy in the reacting codes. Sub- 
routines ANKINR, ACHEMP, ACHEM, MOLIER and PROPR are used 
to describe either an ideal gas initial flow field or an 
equilibrium flow field. Subroutine DIFF is responsible 
for marching- forward the solution as it does for non-reacting 
cases, except the vector U now has more than four components 
in order to consider the chemical species. A subroutine 
DIFFUS is added to account for the mass diffusion vector, 
which is not used in non-reacting codes. Transport co- 
efficients are calculated in subroutine TRANSC , which calls 
either TRANSP for a multicomponent mixture or TRANSK for 
a simple mixture. Thermodynamic properties are calculated 
also in subroutines THERM and JANAF. The subroutine DISSS 
is used for the calculation of the net rate of production 
of species in Eqn (5) . Several other subroutines are 

used to perform the calculation of transport coefficients, 
they are subroutines GAUSS, GJR, LGRNE, SPLNl and SPLN2. 

It should be mentioned that the reacting blunt body codes 
use the same name of subroutines, as in the non-reacting codes, 
for example, MAIN, INIT, OUTPUT, etc., but the content of 
these subroutines are not identical. 
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3.0 PROGRAM USAGE 

3 . 1 Input Description 

3.1.1 Data Specification and Definitions 



NON- 

REACTING THIN 

OR THICK SHOCK CODE 

Card 

Format 

Symbol 

Content 

1 

12A6 

ZA 

Title of the problem 

2 

1415 

NRUN 

MONTH 'I 

Number of the run 



NDAY 1 

MYEAR I 

> Date of the calculation 
1 



) 

LB 

0 for planar flow, 1 for 
axisymmetric flow 



LE 

1 for ellipsoid, 2 for 
paraboloid, 3 for hyperboloid 
(if LB = 1) 

3 

1415 

MA 

Number of mesh spacings along 


the body 

NA Number of mesh spacings across 

the shock layer or the com- 
putational region 

KA Number of time steps to be 

executed 

JA Number of time steps at which 

complete output is required 

4 7E10.4 EPS ( 1) 0.001-0.0001 for the Newtonian 

iteration procedure used in 
SHOCK 
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Card 

Format 

Symbol 

Content 

5 

7E10.4 

RE 

Free stream Reynolds number per foot 



PR 

Free stream Prandtl number 



TW 

Normalized body temperature <30 



BETA 

Parameter for coordinate 




transformation 



RN 

Nose Radius, ft 

6 

7E10.4 

UO 

Free stream speed , ft/sec 



PIN 

2 

Free stream pressure , Ib/ft 



HIN 

2 

Free stream enthalpy, (ft/sec) 



GAMMA 

Free stream ratio of specific 




heats 

7 

7E10.4 

STAB 

0.1-1 stability criterion for 
time increment 



ELL 

Ratio of the major and minor 
axis of an ellipsoid 



ZO 

) Location of the origin relative 



YO 

j to the origin of body coordinates 



XIMAX 

Location of downstream Line 



ANGLE 

Half Angle of Conical Afterbody 


Control cards are needed to delete ACHEMP for an equilibrium 
air calculation or to delete ACHEM for a perfect gas cal- 
culation . 
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REACTING THIN OR THICK SHOCK CODE 


Card 

Format 

Symbol 

Content 

1 

12A6 

2A 

Title of the problem 

2 

1415 

NRUN 

Number of the run 



MONTH ) 

NDAY j 

MYEAR > 

1 Date of the run 
1 



LB 

1 for axisymmetric flow 
0 for planar flow 



LE 

1 for ellipsoid 

2 for paraboloid 

3 for hyperboloid 



JC 

Number of dependent variables 
9 for five species 
11 for seven species 

3 

1415 

MA 

Number of mesh spacings along 
the body 



NA 

Number of mesh spacings across 
the shock layer or the com- 
putational region 



KA 

Number of time steps to be 




executed 



JA 

Number of time steps at which 
complete output is required 

4 

7E10.4 

EPS (1) 

0.001-0.001 used in the 
Newtonian iteration procedure 
in SHOCK 


5 7E10.4 RE 

Pr 
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Card Format Symbol 

LEWIS 

RN 

TW 

BETA 

6 7E10.4 UO 

PIN 

HIN 

GAMMA 

7 7E10.4 STAB 

ELL 

Z° | 

YO I 

XZMAX 

ANGLE 

The following cards are used 
properties 

8 6 (A6 , 6X) DAT ( 1 ) 


DAT ( 2 ) 


Content 

Free stream Lewis number 
Nose radius, ft 

Normalized body temperature <30 
Parameter for coordinate 
transformation 

Free stream speed, ft/sec 

2 

Free stream pressure, lbs/ft 

2 

Free stream enthalpy (ft/sec) 
Free stream ratio of specific 
heats 

Stability criterion for time 
increment 

Ratio of major and minor axis 
of ellipsoid 

Location of the coordinate 
system relative to the origin 
of body coordinates 
Location of downstream line 
Half angle of conical afterbody 

for calculation of transport 

FROZEN for frozen calculation 
of transport properties 
BLANK for equilibrium calcu- 
lation of transport properties 
SIMPLE for simple transport 
coefficients 

BLANK for multicomponent trans- 
port coefficients 


128 


TR2007 


Card 


9 


10 


11-13 


13+4 

x(No. 

of 

Species 

- 1 ) 


Format 


3F10.5 


3A4 , 8X 
3E10.4 


5E15.8 


5E15.8 


Symbol 
DAT (3) 

DAT (4) 

TLOW 

TMID 

THIGH 

NAM ( 1) 
NAM (2) 
NAM (3) 
WM 
FRAC 

HOO 

COEF 
(I , J ,NS) 

COEF 
(I , J ,NS) 


Content 

ISOTHE for isothermal wall 
BLANK for adiabatic wall 
CATALY for fully-catalytic wall 
BLANK for non-catalytic wall 

Defines two temperature 
intervals for the 
temperature coefficients 


Name of the species 

Molecular weight of the species 
Free stream mass fraction of the 
species 

Molar enthalpy at 0°K 

Temperature coefficients for 
two intervals, 1=1,2. Each 
set has J=l, 7 coefficients 
and for species NS. 

Temperature coefficients for 
two intervals, 1=1,2. Each 
set has J=l, 7 coefficients 
and for species NS. 
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Card 

Format 

Symbol 

Content 

13+4 (No. 

5F5.0 

STC (I , J) 

Basic chemical reactions 

of species) 
+ (No. of 



used in the equilibrium 

Reactions) 



transport properties 
calculation 

The following 

cards are 

used for 

calculation of chemical 

rate of production: 



13+4 (No. 

212 

NSPEC 

Number of species 

of species) 
+ (No. of 


N RE ACT 

Number of chemical 

Reactions) 

+1 



reactions 

• 

• 

• 

• 

+ 

to 

12A6 

Z(I) 

Title of the chemical 




reactions formulas 

. . . . + 3 

11 

ITYPE ( I) 

Type of reaction formula 

+NREACT 

8X 




E10.5 

COEFA(I) 



2F10.5 

COEFB (I) 
ETA (I) 

Rate constants 

.... +3 
+N REACT 
+1 




to 

6A6 

PP (IDUM, 

The type of chemical 

... +3+ 
NREACT7 


JDUM) 

reactions 


In addition to the input through reading the cards outlined 
above, these programs have the capability of reading necessary 
parameters from a tape which contains the flow field results 
obtained from a previous computation. 
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3.1.2 Sample Test Data 

Since there is no difference in the input data between the 
thin and thick shock codes, only the non-reacting and reacting 
input data are given separately. 


Non-reacting thin or thick shock code 


IDEAL GAS BLUNT BODY FLOW AROUND A HYPERBOLOID — 10 DEGREES 
0 6 22 72 I 3 

10 LOO 100 50 

0.0C01 0.01 

179282.5 0.71 6. 0. 

20C00. 23.08 2520000. 1.4 

0.5 2.0 0.0 0.0 1.4 10.0 


Reacting thin or thick shock code 


FINITE-RATE FL0WFIELD ABOUT A 2 FT SFPERE ON THE TRAJECTORY POINT NO 3 
9 6 12 72 1 1 ll 

10 15 200 25 

0.0001 0.01 

14400. 0.71 1.5 10. 0. 

25000. 0.04364 1980000. 1.4 

0.7 1.0 0.0 0.0 1.0 10.0 

FROZEN IS0THE CATALY 


131 


1C00. 8000. 30000. 

N 14.008 

.62289649+00 .75996372-03 -.68415018-07 .26029019-11 -.34945540-16 
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3 . 2 Program Run Preparation 
3.2.1 Deck Setup 

Normal control cards for EXEC II of UNIVAC 1108 are required. 
In addition to the JOB, 7/8 FOR , 7/8 XQT and 7/8 EOF cards, 
7/8 ASG is also needed for Unit A. If output to a tape and 
read information from a tape are desired. Unit R and Q are 
required, respectively for such purposes. Some of the 
subroutines have the same name but different compiled name, 
so deletion cards are also used to prevent any duplication 
in the entry point of a subroutine. A source deck setup is 
shown in Figure 1. 
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Figure 1. Deck setup for the viscous blunt body flow field computations 
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3.2.2 Required I/O Devices 

Only standard input/output units are required. For 
example they are units A, or B, Q and R. 

3. 3 Output Description 

All the printed-out data are in the non-dimensional form. 

The nondimensionalization of these data can be restored by 
using the relations in Section 2.2 of Part I. Data are 
available at the final time step, or at any step by using 
the parameter JA on the input card. Data can also be written 
on a tape if so desired. 

The formats for the data print out are self-explanatory. 

For every mesh point the first line gives the pressure, 
temperature, density, velocity components, specific total 
internal energy, Mach number, location of the mesh point 
measured along the shock standoff distance and the percentage 
of the enthalpy dissipation. The second line is for the mass 
concentrations of species and followed by the local molecular 
weight of the mixture. Note only the first line is printed 
out in the non-reacting codes. 

On the last mesh line corresponding to the surface of the 
body, the coefficients of friction and heat transfer are 
printed out. The percentage of the enthalpy dissipation 
is not shown for in a highly-cooled body it is close to 
unity. 

If the detailed flow field results are not required at 
intermediate time steps, the program writes out a few 
essential parameters, viz., the shock standoff distance, 
the stagnation pressure and the shock speeds. 
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4.0 EXECUTION CHARACTERISTICS 

4 . 1 Restrictions 

Theoretically speaking, the only known restriction on the 
mathematical model used for blunt body flow problems is 
that the free stream Reynolds number should be larger than 
one. Since the Navier-Stokes equations are the second-moment 
equations of the second approximation to the Boltzmann equation, 
the stress and diffusion vectors should have derived from a 
higher- moment theory. Furthermore, the surface non-slip 
velocity and non-jump temperature conditions represent a 
crude condition for highly rarefied gas flow. Thus the 
application of the present model to the transitional flow 
regime must be carried out with great caution. For the other 
limiting case of Reynolds number equal to 10 6 or higher, the 
numerical dispersion and dissipation may dominate the physical 
damping in the flow field computation if the mesh spacing used 
is not sufficiently small. However, if the computer storage 
and cost are not the major concern to the user, the numerical 
scheme developed to solve the mathematical model is adequate 
for nearly all engineering problems. The experiences with 
the viscous reacting program are gained entirely from the 
MSC UNIVAC 1108 system. It has 16 3K (octal) storage and is 
considered to be a rather fast machine for scientific com- 
putations. More than 150 mesh points across the shock layer 
and 20 mesh points along the body can be used for the non- 
reacting codes. A relatively large mesh size must be used 
in the reacting codes, since the number of dependent 
variables are higher and more than ten subroutines are needed 
for the calculation of transport properties and chemical 
rates of production. Nevertheless, 60 points across the 
layer and 15 points along the body can be used for viscous 
reacting flow computation. 
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A preliminary understanding of the problem to be solved 

would be highly beneficial to the economic application 

of the programs. For example, it is not always necessary 

to use the maximum allowable mesh points for the problems 

considered. If the nose region is of interest, 10 mesh 

points along the body is enough. If the free stream 

3 

Reynolds number is less than 10 , the shock layer is fully 
viscous, then 20 mesh points across the layer is adequate. 

It should be noted that if the Reynolds number is less 
than 500 and the thick shock code is used, the number of 
mesh points should be in the neighborhood of 40 , in order 
for the shock structure to be accurately defined. 

The transformation of the first computational plane to the 
second plane to achieve a finer space resolution near the 
body should be used with care. The finer space resolution 
at the body is obtained at the expense of a sparser space 
resolution at the shock. This coordinate transformation, 
therefore, creates difficulties in the calculation of the 
shock. It has been found that for the same number of mesh 
points the rate of convergence is slower when the computa- 
tion is made in the second computational plane instead of 
the first. The rate of convergence toward the steady 
solution becomes worse with a larger value of BETA, which 
provides a higher degree of resolution near the body. The 
recommended combinations of the number of mesh points across 
the shock layer and the value of BETA are N=20 , 6=1 ; N=30, 
6=2; or N=45 , 6=3, etc. 

There are limitations on the free stream speed due to the 
assumptions of vibrational equilibrium and of negligible 
radiative energy transfer. The calculation of enthalpy, 
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specific heat, and the Gibbs functions are valid for 
temperatures lower than 44,000°K. Numerical difficulties 
will occur when the temperature within the shock layer is 
higher than the given upper limit. In the applications of the 
the viscous blunt body program to the trajectory point 4, 
as discussed in Part 1, the flow field results are not as 
accurate as the results for other trajectory points at lower 
altitude, because of the thermodynamic properties of air 
currently in use. 

Finally, some remarks on the selection of coordinate systems 
are necessary. The body intrinsic system provides the 
skin friction and heat transfer coefficients as well as other 
flow properties in the shock layer. But this coordinate 
system can only be used for a smooth body with no discontinuity 
in the curvatures . 

If the flow field beyond the sonic line is of interest, both 
the polar and the cylindrical coordinate system can be used, 
although the polar system is more preferable in practical 
applications. The shortcomings of these two coordinate systems 
is that the calculation of the skin friction and heat transfer 
coefficients is a little involved. 

In conclusion, a sufficient large number of mesh points and 
the body intrinsic coordinates should be used for the accurate 
calculation of body properties in a high Reynolds number flow. 
All three types of coordinate systems and a fewer number of 
mesh points can be used if the major purpose of the compu- 
tation is to obtain flow properties in the flow field. 

These programs can be used efficiently to start the down- 
stream supersonic flow field computations. 
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4 . 2 Running Time and Accuracy 

The running time (or the computation cost) in a finite- 
difference technique is directly related to the accuracy 
of the results. The accuracy obtained is proportional 
to the decreased mesh spacing and increased number of 
time steps. The increase in accuracy is high when the 
mesh spacing is first reduced and the number of time 
steps is increased. But, further improvement of the 
numerical results is not discernable with finer and. finer 
mesh spacing and more and more time steps. Two criteria 
have been found to be useful in determining the accuracy 
and the steadiness of the results. One is the stagnation 
heat transfer coefficient or the largest skin friction 
predicted. The accurate calculation of these values 
depends on the adequate mesh spacing adjacent to the wall. 

A coarse mesh spacing often leads to erroneous values 
mainly because of the resultant poor resolution of 
temperature and velocity profiles near the wall. 

Since the existing flow field analyses in the prediction 
of skin heat transfer and friction coefficients are considered 
to be satisfactory, a comparison between the present and the 
existing results would be helpful. The effect of mesh 
spacing on these coefficients is shown in Figure 5. The 
second criterion concerns the steadiness of the solution. 

It is found that when the thin-shock formulation is used 
the shock speed serves as a good indicator of the time- 
asymptotic solution. The shock speed at first changes 
rapidly from time step to time step. As steady state is 
approached, the shock speed decreases until the shock lo- 
cation is finally stationary and the shock speed maintains 
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a very low value. An example of this is shown in Figure 
4. If the thick-shock formulation is in use, flow property 
profiles plotted at several time steps should be used to 
indicate the steadiness of the solution. The steady-state 
solution is reached when the changes in those profiles are 
not observable. Figure 7a provides the transient solutions 
of temperature and the solution obtained at time step K=1000 
is very close to the steady state result. 
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5.0 REFERENCE INFORMATION 

5 . 1 Program Listing 

I. Non-reacting thin-shock code 
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o:.245 5 7* 50 CONTINUE 

0,245 58* C 

0i,250 59* RETURN 
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FLOWFIEID CALCULATION) FOR ORgITER NOSE OF 5 Ft AT TRAJECTORY POINT 5 
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